Pharmacokinetic–pharmacodynamic modeling of maintenance therapy for childhood acute lymphoblastic leukemia

In the treatment of childhood acute lymphoblastic leukemia (ALL), current protocols combine initial high-dose multiagent chemotherapy with prolonged oral therapy with 6-mercaptopurine (6MP) and low-dose methotrexate (MTX) maintenance therapy. Decades of research on ALL treatment have resulted in survival rates of approximately 90%. However, dose-response relationships vary widely between patients and insight into the influencing factors, that would allow for improved personalized treatment management, is insufficient. We use a detailed data set with measurements of thioguanine nucleotides and MTX in red blood cells and absolute neutrophil count (ANC) to develop pharmacokinetic models for 6MP and MTX, as well as a pharmacokinetic–pharmacodynamic (PKPD) model capable of predicting individual ANC levels and thus contributing to the development of personalized treatment strategies. Here, we show that integrating metabolite measurements in red blood cells into the full PKPD model improves results when less data is available, but that model predictions are comparable to those of a fixed pharmacokinetic model when data availability is not limited, providing further evidence of the quality of existing models. With this comprehensive model development leading to dynamics similar to simpler models, we validate the suitability of this model structure and provide a foundation for further exploration of maintenance therapy strategies through simulation and optimization.

Supplementary Table S1 shows the values of the fixed parameters in the MTX PK models, which were calculated in the following way: For the models , the intercompartmental rates and , the elimination rate and the central volume are the mean of the corresponding five values listed in Table 2  are the corresponding volumes in Table 2 in Ogungbenro et al. converted to L/m 2 by using the covariates WEIGHT and BSA for scaling. [24] Parameters not listed here were taken directly from the literature source and not converted in any way.

C Overview and calculation of fixed PK parameters for the 6MP PK model
Supplementary Supplementary Table S2 shows the values of the fixed parameters in the 6MP PK models, which were calculated in the following way: The bioavailability 6 is taken directly from Brunton et al. [46] The absorption rate 6 , the elimination rate 6 and the central volume 6 are calculated by using the means of in Table 1 in Lennard et al.,[30] the stated dose of 75 mg/m 2 and the explicit solution of the ODE system of a 1-CTM-model: Which implies the following relationships: The absorption rate 6 and the elimination rate 6 were converted to 1/day, and the central volume 6 to L/m 2 .

D Overview of model variants
Supplementary Influx into the red blood cells from the central compartment using Michaelis-Menten kinetics, fixed parameters as in in Supplementary Table S1 online.   MTX1-3, MTX4b   , , , Influx into the red blood cells from the central compartment using Michaelis-Menten kinetics, fixed parameters as in in Supplementary Table S1 online, but with dose-dependent bioavailability as in .

MTX1-3, MTX4b
, , Influx into the red blood cells from the peripheral compartment using Michaelis-Menten kinetics, fixed parameters as in in Supplementary Table S1 online.

MTX1-3, MTX4d
, , , Influx into the red blood cells from the peripheral compartment using Michaelis-Menten kinetics, fixed parameters as in in Supplementary Table S1 online, but with dose-dependent bioavailability as in .

MTX1-3, MTX4d
, , Influx into the red blood cells from the central compartment using Michaelis-Menten kinetics, fixed parameters as in in TSupplementary Table S1 online.

MTX1-3, MTX4b
, , , Influx into the red blood cells from the central compartment using Michaelis-Menten kinetics, fixed parameters as in in Supplementary Table S1 online, but with dose-dependent bioavailability as in .

MTX1-3, MTX4b
, , Influx into the red blood cells from the peripheral compartment using Michaelis-Menten kinetics, fixed parameters as in in Supplementary Table S1 online.

MTX1-3, MTX4d
, , , Influx into the red blood cells from the peripheral compartment using Michaelis-Menten kinetics, fixed parameters as in in Supplementary Table S1 online, but with dose-dependent bioavailability as in .

MTX1-3, MTX4d
, , Influx into the red blood cells from the central compartment using linear kinetics, fixed parameters as in in Supplementary  Table S1 online.

MTX1-3, MTX4a
, , , Influx into the red blood cells from the central compartment using linear kinetics, fixed parameters as in in Supplementary  Table S1 online, but with dose-dependent bioavailability as in .

MTX1-3, MTX4a
, , Influx into the red blood cells from the peripheral compartment using linear kinetics, fixed parameters as in in Supplementary Table S1 online.

MTX1-3, MTX4c
, , , Influx into the red blood cells from the peripheral compartment using linear kinetics, fixed parameters as in in Supplementary Table S1 online, but with dose-dependent bioavailability as in .

MTX1-3, MTX4c
, , Influx into the red blood cells from the central compartment using linear kinetics, fixed parameters as in in Supplementary  Table S1 online.

MTX1-3, MTX4a
, , , Influx into the red blood cells from the central compartment using linear kinetics, fixed parameters as in in Supplementary  Table S1 online, but with dose-dependent bioavailability as in .

MTX1-3, MTX4a
, , Influx into the red blood cells from the peripheral compartment MTX1-3, MTX4c using linear kinetics, fixed parameters as in in Supplementary Table S1 online.   , , , Influx into the red blood cells from the peripheral compartment using linear kinetics, fixed parameters as in in Supplementary Table S1 online, but with dose-dependent bioavailability as in .

MTX1-3, MTX4c
, , , Influx into the red blood cells from the central compartment using Michaelis-Menten kinetics, fixed parameters as in in Supplementary Table S1 online. Parameter without interindividual variability.

MTX1-3, MTX4b
, , , , Influx into the red blood cells from the peripheral compartment using linear kinetics, fixed parameters as in in Supplementary Table S1 online, but with dose-dependent bioavailability as in . Purely additive residual error model.

MTX1-3, MTX4c
, , , , Influx into the red blood cells from the peripheral compartment using linear kinetics, fixed parameters as in in Supplementary Table S1 online, but with dose-dependent bioavailability as in . Purely proportional residual error model.

MTX1-3, MTX4c
, , , , Influx into the red blood cells from the peripheral compartment using linear kinetics, fixed parameters as in in Supplementary Table S1 online, but with dose-dependent bioavailability as in . Initial values with the possibility of varying similarly to the residual error model:

MTX1-3, MTX4c
6MP PK model 6 Influx into the red blood cells using Michaelis-Menten kinetics, fixed parameters as in 6 in Supplementary Table S2 online. 6MP1-2, 6MP3b 6 Influx into the red blood cells using Michaelis-Menten kinetics, fixed parameters as in 6 in Supplementary Table S2 online.

6MP1-2, 6MP3a
, 6 Influx into the red blood cells using Michaelis-Menten kinetics, fixed parameters as in 6 in Supplementary Table S2  Influx into the red blood cells using Michaelis-Menten kinetics, fixed parameters as in 6 in Supplementary Table S2 online. Initial values with the possibility of varying similarly to the residual error model:

PKPD model
The . As the model extensions with adjusted bioavailability and basing the influx on the peripheral compartment did not lead to better results than the base variant with fixed bioavailability and the influx based in the central compartment, , , would be preferred if only looking at this structural model. The second part in Supplementary Table S4 online compares the model variants with  linear influx kinetics, leading to   , , , as the model with the lowest objective function, but with a value of -37,359, it is still higher than the one of model has the lowest objective function value, having a closer look at the results of the parameter estimation revealed that the 95%-confidence interval (CI) for the coefficient of variation of of this model covers 0. We therefore estimated the model again with the interindividual variability for removed, resulting in one common parameter for all individual patients. As the third part in Supplementary  Table S4 online  , without any 95%-CI of the interindividual variabilities covering 0.
We then compared the results of the cross-validation for both of the structurally different model variants for all data and for the cross-validation set, showing only slight differences in the RMSEs and MAEs of the individual predictions, but significant changes in the resulting parameters. Most evidently, is almost two orders of magnitude smaller when estimating with the cross-validation data set, with its relative standard error increasing to 84%. Supplementary Figure S7 online depicts the changes in the fixed effect parameters and by comparing all individual patient parameters resulting from estimating using the whole data set and using the cross-validation data set, showing huge differences between both.
Next, we tested the final model

Residual unexplained variability
Additive residual error 0.000019 (8) 0.000021 (9) Proportional residual error 0.023 (7) 0.0088 (24)  led to a similar objective function value, and the same RMSEs and MAEs of the individual predictions as model 6 , without any 95%-CI of the coefficients of variation covering 0.

Errors of the individual predictions
In the next step, we compared the results of the cross-validation of both models , 6 and 6 . The results of the latter one are presented in Supplementary Table S7 online, comparing the parameter values and the RMSEs and MAEs of the individual predictions of both the estimation using the whole data set and with the crossvalidation set, showing again an expected increase in the RMSEs and MAEs, and almost no changes in the resulting parameters. In contrast to the MTX PK model with linear kinetics, this model leads to significant pvalues < 0.05, rejecting the hypothesis that the true mean of the -estimates is 0. Supplementary Figure S13 online plots the individual patient parameters of model 6 of the estimation using the whole data set against the cross-validation, showing only minor differences between the two.

Errors of the individual predictions
In contrast to the MTX PK models, the trajectories resulting from model We then tested the model ). The model variant with the initial value estimated again led to a result with implausible parameter values and NONMEM not being able to estimate the covariance matrix and is thus not included here. All three other models resulted in a higher objective function value than

G Detailed Results for the full PKPD model
Supplementary To be able to compare the objective function values of the different models for the model evaluation, we estimated the parameters for both PK submodels for all PKPD models, irrespective of them both being included in the effect function. The RMSEs and MAEs of the individual predictions show only negligible differences, with the RMSEs and MAEs of the PK submodels not differing from the results when estimating only the corresponding PK models.
The objective function values are comparable for both models with only linear kinetics in the PK submodels ( 6 and 6 , ), and for both models with linear kinetics in the MTX PK model and Michaelis-Menten kinetics in the 6MP PK submodel ( ), hence including E-MTX in the effect function does not lead to better results than basing the effect function only on E-TGN. Overall, the models with Michaelis-Menten kinetics in the 6MP PK submodel led to lower objective function values than the two models with purely linear PK kinetics. Comparing all models with regard to their ability to capture the full range of observations (with a median of patientwise minima of 0.5 G/L and a median of patientwise maxima of 4.4 G/L) reveals similar results: all models are able to predict the median ANC values well, but display a smaller range for the individual predictions than the observations show, leading to higher patientwise minima (with a median of patientwise minima of 1.14-1.23 G/L) and lower patientwise maxima (with a median of patient-wise maxima of 2.58-2.78 G/L, cf. Supplementary Table S8 online).
The cross-validation of the model 6 led to an immense increase in RMSEs and MAEs of the individual predictions, with mean( ) = 13720 and mean( ) = 3804. Having a closer look at the individual patients revealed an implausible increase in ANC values for one patient with an estimated individual parameter , 6 = 1.02 which was more than three times the size of the population parameter 6 = 0.27, but still in agreement with the model results. The treatment schedule of this patient showed that in the cross-validation data set, the highest 6MP dose of this patient was not included and led to a strong increase in E-TGN values and therefore a severe decrease in ANC values, which then induced the sharp increase in ANC values due to the feedback mechanism (cf. Supplementary Figure S22 online). We therefore concluded that for 6MP, the PK model with linear kinetics is not able to capture the dynamics of 6MP treatment well enough for all plausible treatment schedules and excluded both 6 and 6 , model from further analysis. , the fixed effect parameter 6 still decreases, but again not as sharply as when estimating only the 6MP PK model. The same applies to its relative standard error, which does not show the immense increase for the cross-validation of the full PKPD model. The coefficients of variation show only negligible changes for the cross-validation. Compared to model , 6

Supplementary
, the parameter estimation of the full PKPD model leads to similar results as those of the estimation of only the PK submodels (cf. Tables 2 and 4), except for the 6MP PK submodel, where again both 6 and 6 increase distinctly. Again, integrating E-MTX into the effect function does not improve the model. Figure S20: Goodness-of-fit plots of observed and by the model , 6

Supplementary
predicted E-MTX, E-TGN and ANC values, estimation using the whole data set.

J Exemplary individual patient trajectories
Supplementary Figure S23: Individual ANC trajectories of four patients predicted by the final PKPD model.
Four patients exhibiting qualitatively different oscillating behavior were selected after a visual check of patient ANC trajectories.

K Sensitivity Analysis
Supplementary